Creation and Reproduction of Model Cells with Semipermeable Membrane 



Hidetsugu Sakaguchi 
Department of Applied Science for Electronics and Materials, 
Interdisciplinary Graduate School of Engineering Sciences, 
Kyushu University, Kasuga, Fukuoka 816-8580, Japan 

A high activity of reactions can be confined in a model cell with a semipermeable membrane 
in the Schlogl model. It is interpreted as a model of primitive metabolism in a cell. We study 
two generalized models to understand the creation of primitive cell systems conceptually from the 
view point of the nonlinear-nonequilibrium physics. In the first model, a single-cell system with 
a highly active state confined by a semipermeable membrane is spontaneously created from an 
inactive homogeneous state by a stochastic jump process. In the second model, many cell structures 
00 ' are reproduced from a single cell, and a multicellular system is created. 
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I. INTRODUCTION 



o 

' Reaction-diffusion systems are important models for various nonlinear processes in biological systems. [l|, H, Q Most 
reaction-diffusion systems have been studied in homogeneous media to simplify the processes and their analyses. 
Recently, reaction-diffusion systems in more complicated media such as membranes and microemulsions have been 
/**s \ investigated. [1, [1] 

■ Complicated biochemical reactions take place in a cell. Several researchers considered that chemical reactions in a 
confined cell- like structure are an important step to life. Oparin considered that "coacervate" played an important 
role in the origin of a cell in prebiotic chemical evolution. [6[ He pointed out the importance of chemical reactions 
confined in the cellular structure "coacervate". Dyson pointed out the importance of the jump to a highly activated 
state in a certain bistable system as the origin of life. [3] Recently, Noireaux and Libchaber have studied a cell-like 
bioreactor using a phospholipid vesicle. [f| 

I , A cell is separated from the outside world by a membrane. A cell becomes an independent element owing to 

■ the existence of the membrane. Complicated chemical reactions occur only inside of the cell. Some materials and 
reactants are transported through the membrane. The consumption of nutrients and elimination of waste matter 
through the membrane is a basic process of metabolism. A highly active nonequilibrium state is maintained only 
inside of the cell. Inside and outside of the cell, materials are in an aqueous solution, but the membrane is composed 
of lipid. Therefore, the transport coefficient in the membrane is different from that inside and outside of the cell. We 
constructed a simple reaction-diffusion system based on the Schlogl model to understand the confined nonequilibrium 
states accompanying the transport of materials through the membrane qualitatively. @ In the model, the Schlogl 
reaction proceeds in a vesicle with a semipermeable membrane, where the diffusion constant of some materials is small 
or zero. The semipermeable membrane is necessary to confine the activity inside of the cell. Because complicated 
reactions including those of DNA and RNA, and active transport through the membrane were not included in the 
model, our model is too simple as a model of a living cell. Our model is not a realistic model, but an artificial model 
for considering "metabolism" conceptually from the view point of a dissipative structure far from equilibrium. In this 
paper, we propose two models for generalizing the previous model. In the first model, a single-cell system with a highly 
active state confined by a semipermeable membrane is created from an inactive homogeneous state by a stochastic 
jump. In the second model, semipermeable membranes are reproduced and a multicellular system is created. 



II. SCHLOGL MODEL 



The Schlogl model is a set of chemical reactions including three kinds of chemicals, namely, U, V and W [T(| 
represented as 

V + 2U # 3U, U # W. 
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We interpret the three chemicals V, U and W as a nutrient material, a self- catalytic product and a waste material, 
respectively. The reaction-diffusion equations for u = [U] , v = [V] , and w — [W] are expressed as 

du 

— = k\vu 2 — A^it 3 — k^u + k^w + D u S7 2 u, 
at 

dv 01 1 

— = -k 1 vu 2 + k 2 u i + D v V 2 v, 
at 

duo 

= k^u — k^w + D w \7 2 w. (1) 

where kx, k 2 , &3 and k^ are the rate constants for the reactions V+2U— >3U, 3U— >V+2U, U— >W, and U, 
respectively, and D u , D v , and D w are the diffusion constants of the materials. In general reaction-diffusion equations, 
the diffusion constants are assumed to be uniform in space. However, we assume a vesicle surrounded by a membrane 
as a model of a primitive cell. The diffusion constants inside of the vesicle are different from those in the bulk region. 
We assume nonuniform diffusion constants as 

D u (r) = D u i, for r < r , r > r + d, 

= D uQ for r < r < r + d, 

D v (r) = D vl , for r < r , r > r + d, 

= D v o for r < r < r + d, 

D w (r) = D wl , for r < r , r > r + d, 

= D w0 for r < r < r + d, (2) 

where r is the distance from the center of a spherical cell. The value ro + d represents the size of the cell, and d is the 
thickness of the membrane. The diffusion constants are uniform inside and outside of the cell. The diffusion constants 
inside of the cell might be different from those outside of the cell, but we assume the same values for the sake of 
simplicity. In all simulations in this paper, we have assumed D u \ = D v x = D w \ = D\ and D u q = D v q = D w q = Dq. 
The diffusion constants are assumed to be smaller in the membrane region tq < r < r$ + d. That is, the transport 
of materials by diffusion is assumed to be more difficult in the membrane region. In our previous research, we have 
studied several other cases. In the case of D u q = 0, D v q = D w q =/= 0, a confined active state was obtained, which is 
similar to that in the case in this paper. In another case of D v q = 0, D u q = D w q =/= 0, no confined active state was 
obtained because the nutrient cannot be supplied to the cell. 

Figure 1(a) shows an example of a confined active state at k\ = k 2 = k% = k^ = l,Di = 1,Dq = 0.01, ro = 2, 
and d = 1. The boundary conditions are assumed to be u(R) = w(R) = and v(R) = vq = 0.8 at a radius R = 20 
located far from the center of the cell. The initial conditions are v(r) = Vo,w{r) = 0, and u(r) — for r > ?*o + d, 
and u(r) = 0.07 for r < tq + d. By assuming a spherical symmetry, V 2 is replaced by d 2 jdr 2 + (2/r)(d/dr) in the 
numerical simulation. A high activity of the self-catalytic product U is maintained inside of the cell owing to the 
low diffusivity of the membrane. @ Figure 1(b) shows another stationary state called the death state: v(r) = Vq and 
u{r) = w(r) = 0, because no self-catalytic product is produced and the concentrations u, v and w take the same values 
inside and outside of the cell membrane. This state was obtained from the initial conditions: v(r) = Vo,w(r) = 0, 
and u(r) = for r > ro + d and u(r) = 0.06 for r < ro + d. The active and death states are bistable when Dq = 0.01 
and Dx = 1. When Dq = Dx = 1, only the death state appears if the other parameters are fixed. 

In a living cell, a material Z such as lipid, from which the cell membrane is constructed, is synthesized by the cell 
itself. We can construct a simple model in which the material generating the low diffusivity of the cell membrane is 
synthesized via the self-catalytic product U as 



dz c 
~dt = 47r/3(r + df J 



ro+d 

9(u - u t h)^KT 2 dr - k z z, (3) 



where z is the concentration [Z] of Z, uth is the threshold of u for the creation of the material Z, and 6(x) is the 
Heaviside step function. Here, we assume that Z is produced when u is larger than uth with a constant rate c and is 
decomposed with a rate constant k z just as in a simple model. The diffusion constant D\ = 1 and Dq is assumed to 
linearly decrease as Do = 1 — 0.99z with the existence of the material Z such as lipid. This is also a simple model. 
We have performed numerical simulation for c = l,u t h = 0.03, and k z = 1 from the initial condition: z = 0, Do = I, 
v(r) — vq = 0.8, w(r) — 0, and u(r) = for r > ro + d and u(r) = 1 for r < ro + d. That is, we have assumed a high 
activity of U inside of the cell as an initial condition. Figure 2(a) shows the time evolution of Do- As the material 
Z is produced by a high activity of U inside of the cell, Do decreases from 1 towards 0.01. The low diffusivity of the 
membrane maintains the high activity of U inside of the cell. The low diffusivity of the membrane is maintained by 



(a) 



(b) 




1 2 3 4 5 6 7 




FIG. 1: (a) Profiles of u(r),v(r) and w(r) for the active state obtained using eq. (1) for ki — k^ — k$ — k± — l,Di — I, Do = 
0.01, xo = 2, and d = 1. The initial conditions are v(r) — 0.8, w(r) = 0, and u(r) = for r > ro + d and u(r) — 0.07 for 
r < ro + d. (b) Death state obtained from the initial conditions: v(r) = 0.8, w(r) = and u(r) = for r > ro + d and 
u(r) = 0.06 for r < ro + d. 
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FIG. 2: (a) Time evolution of the diffusion constant Da(t) of the membrane, (b) Localized solutions: u(r), v(r), and w(r) to 
eqs. (1) and (3) with u t h = 0.03. 



the high density of the self-catalytic product U. If Dq is close to 1 for a long time, the concentration of U inside of 
the cell decreases to 0, which leads to the death state. In the death state, z = and Dq remains to be 1. Figure 2(b) 
shows an active state that appeared as a result of the numerical simulation. The profiles of u, v and w are almost 
the same as those in Fig. 1(a). In this model, the low diffusivity of the membrane and the high activity of U inside of 
the cell are mutually dependent on each other. In other words, an active state confined in the membrane with a low 
diffusivity does not appear naturally from a death state in a system of uniform diffusivity Dq = D\. 



III. STOCHASTIC SCHLOGL MODEL 



For the comparison with the stochastic model introduced later, we first consider a one-dimensional model where V 2 
in eq. (1) is replaced by d 2 /dx 2 . Figures 3(a) and 3(b) show an active state and a death state in the one-dimensional 
model for k\ = tt2 = 8 x 10~ 6 ,/s3 = fej = 2, D\ = 2, Dq = 0.02, xq = 2, and d = 1. The boundary conditions are 
du/dx = 0,v = 400, and w = 1 at x = 10, and du/dx = dv/dx — dw/dx — at x — 0. The two states are bistable 
similarly to that in the case of the spherically symmetric model shown in Fig. 1. The active state confined by the 
membrane of low diffusivity is also not naturally created from a death state in this deterministic model. 

We construct a stochastic version of the Schlogl model as a model of the birth of a primitive cell. Stochastic 
chemical reactions have been studied as nonlinear-nonequilibrium systems far from equilibrium. [l| The stochasticity 
originates from the finiteness of the number of molecules. A novel state generated by the discreteness of molecules 
was reported in a small autocatalytic system. [ll| Recently, several stochastic models have been studied to understand 
fluctuations in gene networks. [l2l [l3l| 

The one-dimensional space is further simplified to a discrete one-dimensional lattice of interval 1, and the time is 
also discretized. The numbers of the molecules U, V, and W at the ith site and the nth step are respectively denoted 
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FIG. 3: Bistable states in a one-dimensional model, (a) Active state, (b) Death state. 
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FIG. 4: (a) Time evolution of N u (0). (b) Time averages u,v, and w of N u (i), N v (i), and N w (i). (c) Probability distributions 
P U {N),P V (N) and P W {N) of N u (i),N v (i) and N w (i) at i = 1. (d) Probability distributions P„(N),P V {N), and P m (iV) of 
N u (i),N v (i), and N m (i) at i = 7. P U (N) and P m (N) almost overlap in (c) and (d) . 



by iV u (z, n), A^(i, n), and N w (i, n) where i = 0, 1, 2, • • ■ , M. The probability for the reaction V+2U— >3U to take place 
at the ith site and nth step is assumed to be kiN v (i,n)N 2 (i,n). That is, the numbers of U and V change by 1 as 
N u (i,n + 1) = N u (i,n) + 1 and N v (i,n + 1) = N v (i,n) — 1 with the probability kiN v (i,n)N 2 (i,n). The numbers 
of U and V do not change with the probability 1 — kiN v (i,n)N 2 (i,n). Similarly, the probability for the reaction 
3U^V+2U is assumed to k2N u (i,n) 3 , and the probabilities of the reaction U — > W and the inverse reaction are 
respectively k^N u {i,n) and ki,N w {i,n). We perform a Monte-Carlo simulation in which the numbers of molecules, 
namely N Ul N v , and N w exhibit random walks according to the probabilities determined by the molecule numbers at 
each site. The diffusion processes are also simulated by another random walk between neighboring sites. That is, the 
number N u (i) of U at the ith site increases by 1 as N u (i, n + 1) = N u (i, n) + 1, and N u (i + 1) at the (i + l)th site 
decreases by 1 as N u (i + l,n + 1) = N u (i + 1, n) — 1 with the probability D u (i + \/2)N u (i + 1, n). Inversely, N u (i) 
decreases by 1 as N u (i, n + 1) = N u (i, n) — 1 and N u (i + 1) increases by 1 as N u (i + 1, n + 1) = N u (i + 1, n) + 1 
with the probability of D u (i + l/2)N u (i,n). Each diffusion constant D u (i + 1/2) is defined at the bond between 
the ith site and the (i + l)th site. The diffusion constant is assumed to be D u (i + 1/2) = D\ except for i = 2 and 
D u (i + 1/2) = Dq only for i = 2. At the left edge i — 0, D u (i + 1/2) is assumed to be 0. The reaction and diffusion 
processes are performed independently in our Monte-Carlo simulation. 

If the fluctuations of N U ,N V , and N w are neglected, that is, N U ,N V , and N w are assumed to be equal to their 
ensemble averages (N u ), (N v ) and (N w ), similarly, N^N V and are assumed to be equal to (N U ) 2 {N V ) and (N u ) 3 , 
rate equations for iV u , N v , and N w are obtained as 

(JV u (i,»+l)) = (N u (i,n)) + k 1 {N v (i,n)){N u (i,n)) 2 - k 2 {N u (i,n)} 3 - k 3 (N u (i,n)) +k i (N w {i 1 n)) 
+ D u (i+ l/2,n)((N u (i + l,n)) - (N u (i,n))) + D u (i - l/2,n)((JV u (t - l,n)> - (JV u (i,n))), 

(JV„(i s n + l)) = (iV v (i, n)) — ki{N v (i, n)){N u (i, n)} 2 + k 2 {N u {i, n)} 3 
+ D v (i + 1/2, n){(N v (i + 1, n)> - (iV„(i, n)» + D„(i - 1/2, n)((iV w (i - 1, n)> - (^(i, »))), 

(N w (i,n + 1)) = (N w (i,n)) + k 3 {N u (i,n)) - k A {N w (i,n)) 
+ D w (i + l/2,n)({N w (i + l,n)) - (N w (i,n))) + D w (i - l/2,n)((N w (i - l,n)> - (iV w (i,n))). (4) 
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FIG. 5: (a) Time evolution of N u (0). (b) Time evolution of Z. (c) Long-time averages of N u (i), N v (i), and N w (i). 



If the difference operators in both time and space in eq. (4) are replaced with a differential operator, the original rate 
equation eq. (1) is recovered. 

We have performed a direct Monte-Carlo simulation for fci = k-i = 8 x 10 -10 , k% = k± — 2 x 10~ 4 , D\ = 2 x 10~ 4 , and 
Dq = 0.02 x 10~ 4 , and N v and N w are respectively fixed to be N v = 400 and N w = 1 at the boundary site i = M = 10. 
The relative ratios of these parameters are the same as those in the case of the deterministic one-dimensional model 
shown in Fig. 3. The initial condition was a death state, that is, N u (i) = N w (i) = and N v (i) = 400. Figure 4(a) 
shows the time evolution of N u (0) at the left edge site. N u (0) increases rapidly at t = 1300. Here, the time t is 
determined using t — nAt, where At = 10~ 4 . It is a stochastic jump process from a death state to an active state. We 
observed no inverse transition from the active state to the death state until the final time of the simulation t = 10000. 
Figure 4(b) shows long-time averages of N u (i), N v (i), and N w (i) after a transition time. The profiles are fairly close 
to those for the deterministic model shown in Fig. 3, although the space is continuous in the deterministic model. 
Figures 4(c) and (d) show the probability distributions P U (N), P W (N), and P V (N) for N u (i), N w (i), and N v (i) at 
two sites i = 1 and i — 7. The probability distribution of P W (N) almost overlaps with P U (N) in these figures. The 
stochastic transition has occurred owing to the fluctuations in N u , N v , and N w . 

We can further construct a model in which the material Z producing the low diffusivity of the membrane is 
synthesized by the self-catalytic product U as 

2 

Z{n + 1) = Z{n) + At{6{J2 N u (i, n) ~ u th ) - Z(n)}, (5) 

where At = 10~ 4 and u t h is the threshold of u for generating the material. The diffusion constant Dq in the membrane 
region is assumed to be given by Do = Z?i(l — 0.99Z(n)). The quantity Z fluctuates in time because of the temporal 
fluctuation in N u . The initial condition is Dq = D\ = 2 X 10~ 4 , and therefore Z(0) — 0. At the boundary site 
i = M = 10, N v and N w are respectively fixed to be N v — 400 and N w — 1. The initial condition is a death 
state N u (i) = N w (i) — and N v (i) — 400. (Almost the same result was obtained from another initial condition 
N u (i) = N w (i) = N v (i) = 0.) Figures 5(a) and 5(b) show the time evolution of A^ u (0) at the left edge i = and Z(n) 
for uth = 8. Initially, N u (0) exhibits fluctuation around 1, and Z increases intermittently from 0, but it decreases to 
nearly again for t < 27500. However, N u (0) increases rapidly to N u = 130 at t ~ 27500, and then Z increases to 1. 
After that, Z is fixed to 1 and then Dq is fixed to 0.01-Di. Figure 5(c) shows long-time averages of N u (i), N v (i), and 
N w (i) after a transition time, which is almost equal to those in Fig. 4(b). In a uniform system with Dq = D\, only 
the death state is stable. Therefore, the semipermeable membrane with low diffusivity and the active state inside 
of the cell should be created almost simultaneously. Actually, such an active state surrounded by a semipermeable 
membrane is spontaneously created from a death state in a system with an almost homogeneous diffusivity by a 
stochastic jump process at t ~ 27500. This stochastic jump is more difficult than the simple stochastic jump in a 
bistable system studied in the previous model shown in Fig. 4. The high activity inside of the cell is maintained and 
protected by the semipermeable membrane, and the membrane that determines the form of the cell is created by the 
materials inside of the cell. Both are necessary for the construction of a cell. Therefore, the stochastic jump process 
might be interpreted as the birth of a cell, and the jump process might also be interpreted as one of the sequential 
jump processes toward the origin of life, as Dyson has pointed out. 
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FIG. 6: Snapshot profiles of D(x, y) at (a) t — 5 and (b) t = 150 for vq = 1.2. The shaded regions represent the semi-permeable 
membrane where the diffusion constant D(x,y) is smaller than 0.4. (c) Time evolution of D(x,y,t) at the section of y = L/2. 




FIG. 7: Snapshot profile of D(x,y) at (a) t = 5 and (b) t — 75 for wo = 0.7. The shaded regions represents the semi-permeable 
membrane where the diffusion constant D(x,y) is smaller than 0.4. (c) Time evolution of D(x,y,t) at the section of y = L/2. 



IV. REPRODUCTION AND EXTINCTION OF MODEL CELLS 



In §2, we have constructed a single cell model with a spherical symmetry. The low diffusivity in the membrane 
region is maintained by the material Z produced inside of the cell. We can construct a two-dimensional model where 
V 2 in eq. (1) represents d 2 /dx 2 + d 2 /dy 2 . Each cell is assumed to have a square form of size I = 2, and the width d 
of the membrane for each cell is assumed to be 0.2. The diffusion constant Dq in the membrane region for each cell 
at site is given by Do(i,j) = (1 — 0.99z(i, j))D\. The material Z at the site is assumed to be produced by 
the cell surrounded by the membrane as 

dz(i,j) 1 f u+1)l 

— ^ — = p J J 9(u(x,y) -uth)dxdy - z(i,j). (6) 

We performed numerical simulations of eqs. (1) and (6). The parameters are k\ = k,2 = k% = k± = 1 and D± — 1. The 
boundary conditions are v = vq,u = w = for r > 21. The nutrient concentration vo at the boundary is a control 
parameter. The initial condition is z(i, j) = 0, v(x, y) — Vo, w(x, y) — 0, and u(x, y) = 1 for r < lo and u(x, y) — for 
r > lo, where r = ^(x - L/2) 2 + (y - L/2) 2 < 2 is the distance from the center of the system of size L. That is, only 
the central region of size Iq is activated as an initial condition. The membrane is created at the central site by the 
existence of U. If the membrane is constructed around the cell region, the concentration of the self-catalytic product is 
maintained inside of the membrane, and a single cell is created at the central site. By the slow diffusion of U through 
the membrane, U increases in neighboring sites, leading to the construction of the membrane in neighboring sites, 
and new cells are created around the central cell. In this deterministic model, an active cell cannot be created from 
a death state as in the previous section, but a multicellular structure can appear starting from a single active cell. 

Figures 6(a) and 6(b) show two snapshots of cellular structures at t = 5 and 150 for vq = 1.2, Iq = 2, and L = 42. 
The membrane region with a small diffusion constant Do ~ O.OlDi is marked in Fig. 6(a). Nine cells are created at 
t = 5, and the cell number increases with time, as shown in Fig. 6(b). Figure 6(c) shows the time evolution of D(x, y) 
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at the section y = L/2. The cell number increases monotonically. Similar self-replicating structures were also studied 
in several models of artificial life. [3] Figures 7(a) and 7(b) show two snapshots of the cellular state at t = 5 and 75 
for a smaller nutrient concentration vq — 0.7 and Iq — 2.5. The system size is L = 30. Figure 7(c) displays the time 
evolution of D(x,y) at the section y — L/2. Because the initially activated region is slightly larger than that in the 
case of Fig. 6, several cells are created at t = 5. However, the cell number decreases with time and finally the cell 
structure disappears, which leads to a death state: u = w = for vq = 0.7. The extinction is caused by the lack of 
the nutrient V at the outer boundary. 



V. SUMMARY 



The membrane works as the boundary from the outside world, which makes a cell an independent element separated 
from the outside world. The membrane as a container is constructed by the contents of the cell, and a high activity 
of the contents of the cell is maintained and protected by the membrane. We have constructed a stochastic Schlogl 
model to understand the creation of the mutually dependent relation. We have shown with a Monte-Carlo simulation 
of the model that an active state surrounded by a semipermeable membrane is self-organized via a stochastic jump 
process starting from a death state in a homogeneous system with respect to diffusivity. It can be interpreted as 
the birth of a cell. We have also constructed a breeding model of cell structures. If the nutrient is sufficiently large, 
the cell number increases with the creation of the membrane, which leads to a multicellular system. If the nutrient 
is not sufficient, the cell structure is extinguished. Our model is a toy model, but might be a suggestive model for 
considering a primitive cell or life from the view point of dissipative structures far from equilibrium. 
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